﻿(* PolyCopm.fs includes the implementation of Polynomial Complexity algorithm.
   The algorithm pseudo-code can be found in "Polynomial Complexity.docx"       *)
   

module POLY_COMP_ANALYSIS
open Microsoft.FSharp.Math
open Microsoft.FSharp.Collections
open BJK
open FSEvent


(* ------     -----------    --------- --------- --------- --------- *)
(* ------     -----------    --------- --------- --------- --------- *)
(* ------     Abstract semantics definitions     --------- --------- *)
(* ------     -----------    --------- --------- --------- --------- *)
(* ------     -----------    --------- --------- --------- --------- *)


(* D it the basic type used by Polynomial Complexity algorithm.      *)
(* D represent variable degree relation.                             *)
type D = Degree of int | Epsilon

(* An Abstract Monomial Bound(AMB) is a D vector where AMB[i]        *)
(* represent the degree of Xi                                        *)
type AMB = D list

(* An Abstract Polynomial Bound(APB) is a set of AMBs.               *)
type APB = AMB list

(* BPR stands for Basic Polynomial Relation and defines a relation   *)
(* of a command's output variables to its input variables,           *)
(* BPR is a n*n matrix made of n column vectors for AMBs'            *)
(* representing variable Xi, 0<i≤n.                                  *)
type BPR = Matrix<D>

(* PR stands for Polynomial Relation and defines a set of BPRs       *)
type PR  = BPR list 



(* ------     -----------    --------- --------- --------- --------- *)
(* Generic Max                                                       *)
(* Max(List, >= operator) : subset List                              *)
(* ------     -----------    --------- --------- --------- --------- *)
(* The Max function inputs are a list "listInput" of elements and an
   order operator.
   The Max function returns a new list "listOutput" where
   listOutput⊆listInput, for every x∈listInput exist y∈listOutput
   such that x≤y and no strict subset of listOutput satisfied same
   requirements.                                                     *)
(* ------     -----------    --------- --------- --------- --------- *)
(* Notice that two elements in the List might not ordered,           *)
(* Example: Max{[1,1], [1,2], [2,1]} = {[1,2], [2,1]}                *)
(*   because [1,2]>=[1,1] (and [2,1]>=[1,1])                         *)
(*   but [2,1]!>=[1,2] and [1,2]!>=[2,1]                             *)
(* ------     -----------    --------- --------- --------- --------- *)
let rec Max(listOfElem : 'a list, greaterEqualFunc) : 'a list =
    (* RemoveSmaller return the subset - list minus {elem | elem of list and elemToCheck>=elem} *)
    let inline RemoveSmaller(listOfElem : 'a list, elemToCompare : 'a, greaterEqualFunc) = snd (List.partition (fun elem -> greaterEqualFunc elemToCompare elem) listOfElem)
    match listOfElem with
        | []                          -> []
        | (first_elem::[])            -> [first_elem]   // Max of list with one element is the list itself  
        | (first_elem::rest_of_list)  -> 
            if List.exists (fun elem -> greaterEqualFunc elem first_elem) rest_of_list                // Is there an element in the rest of the list that is "greater or equal" from first_elem
                then Max(rest_of_list, greaterEqualFunc)                                                         // No need for first_elem, Use the rest of the list
                else first_elem::Max(RemoveSmaller(rest_of_list, first_elem, greaterEqualFunc), greaterEqualFunc) // Save first_elem and append with subset created form {Rest of the list}/{elements that are smaller of equal to first_elem}


(* ------     -----------    --------- --------- --------- --------- *)
(* Type D operation definitions                                      *)
(*                                                                   *)
(* We will use the character ? to mark D operators                   *) 
(* ------     -----------    --------- --------- --------- --------- *)

let rec ( *? ) (d1 : D) (d2 : D) : D = 
    match (d1, d2) with 
        | (Epsilon      , Degree(0) )   -> Degree(0)
        | (Epsilon      , _         )   -> Epsilon
        | ( _           , Epsilon   )   -> d2 *? d1

        | (Degree(d1)   , Degree(d2))   -> Degree(d1*d2)

let rec ( +? ) (d1 : D) (d2 : D) : D = 
    match (d1, d2) with 
        | (Epsilon      , Degree(0) )   -> Epsilon
        | (Epsilon      , Epsilon   )   -> Epsilon
        | (Epsilon      , _         )   -> d2
        | ( _           , Epsilon   )   -> d2 +? d1

        | (Degree(d1)   , Degree(d2))   -> Degree(d1+d2)

(* for easier Order Relation we convert D to float *)
let DtoFloat(d : D) : float = 
    match d with 
        | (Degree(d)  )   -> float d
        | (Epsilon    )   -> 0.1

let inline ( <=? ) (d1 : D) (d2 : D) : bool = DtoFloat(d1) <= DtoFloat(d2)
let inline ( ==? ) (d1 : D) (d2 : D) : bool = d1.Equals(d2)
let inline ( >?  ) (d1 : D) (d2 : D) : bool = not (d1 <=? d2)
let inline ( >=? ) (d1 : D) (d2 : D) : bool = d2 <=? d1


let DtoStr (d : D) : string =
    match d with
        | Epsilon   -> "e"
        | Degree(d) -> string d

        

(* ------     -----------    --------- --------- --------- --------- *)
(* AMB operation definitions                                         *)
(*                                                                   *)
(* We will use the character & to mark AMB operators                 *)
(* ------     -----------    --------- --------- --------- --------- *)
 
let inline ( *& ) (d_mul : D)  (amb  : AMB) : AMB  = List.map (fun d -> d_mul *? d) amb
let inline (==&)  (amb1 : AMB) (amb2 : AMB) : bool = List.forall2 (fun elem1 elem2 -> elem1 ==? elem2) amb1 amb2
let inline (!=&)  (amb1 : AMB) (amb2 : AMB) : bool = not (amb1 ==& amb2)
let inline (<=&)  (amb1 : AMB) (amb2 : AMB) : bool = List.forall2 (fun elem1 elem2 -> elem1 <=? elem2) amb1 amb2
let inline (>=&)  (amb1 : AMB) (amb2 : AMB) : bool = List.forall2 (fun elem1 elem2 -> elem1 >=? elem2) amb1 amb2
let inline (<&)   (amb1 : AMB) (amb2 : AMB) : bool = (amb1 !=& amb2 && amb1 <=& amb2);
let inline (+&)   (amb1 : AMB) (amb2 : AMB) : AMB  = List.map2 (fun d1 d2 -> d1 +? d2) amb1 amb2

(* Set D value at AMB *)
let inline ambSet(amb : AMB, pos : int, d_set : D) = List.mapi (fun i d -> if (i.Equals pos) then d_set else d) amb


(* ------     -----------    --------- --------- --------- --------- *)
(* APB operation definitions                                         *)
(*                                                                   *)
(* We will use the character % to mark APB operators                 *)
(* ------     -----------    --------- --------- --------- --------- *)

let inline MakeAPBfromAMB(arr : AMB) : APB = [arr];

let inline ( *% ) (apb1 : APB) (apb2 : APB) : APB = List.map2 (fun amb1 amb2 -> amb1 +& amb2) apb1 apb2

(* Max of join *)
let inline (|=|) (apb1 : APB) (apb2 : APB) : APB = 
    Max(apb1 @ apb2, (>=&))

    

(* ------     -----------    --------- --------- --------- --------- *)
(* BPR operation definitions                                         *)
(*                                                                   *)
(* We will use the character ^ to mark BPR operators                 *)
(* ------     -----------    --------- --------- --------- --------- *)

(* *^ is D matrix multiplication *)
let ( *^ ) (bpr1 : BPR) (bpr2 : BPR) : BPR = 
    if (bpr1.NumRows <> bpr1.NumCols || bpr2.NumRows <> bpr2.NumCols || bpr1.NumRows <> bpr2.NumRows) then Util.fail("PolyCompAnalysis.BPR*: wrong dim?.")
    else
        let dim = bpr1.NumRows
        let mutable bpr = Matrix.Generic.create<D> dim dim (Degree(0))
        for j in 0 .. dim-1 do
            for i in 0 .. dim-1 do
                for k in 0 .. dim-1 do
                    bpr.[j,i] <- bpr.[j,i] +? (bpr1.[k,i] *? bpr2.[j,k])
        bpr


let inline (==^) (bpr1 : BPR) (bpr2 : BPR) : bool = Matrix.Generic.foralli (fun i j d1 -> d1 ==? bpr2.[i,j]) bpr1
let inline (<=^) (bpr1 : BPR) (bpr2 : BPR) : bool = Matrix.Generic.foralli (fun i j d1 -> d1 <=? bpr2.[i,j]) bpr1
let inline (>=^) (bpr1 : BPR) (bpr2 : BPR) : bool = Matrix.Generic.foralli (fun i j d1 -> d1 >=? bpr2.[i,j]) bpr1

(* We'll create the identity matrix our self since Matrix.Generic.identity is not working very well. *)
let IdentityMatrix(numVars : int) : BPR = Matrix.Generic.mapi<D> (fun i j d -> if (i.Equals(j)) then Degree(1) else d) (Matrix.Generic.create<D> numVars numVars (Degree(0)))

let rec BPRtoSTR(bpr : BPR) : string = 
    let mutable res = ""
    let dim = bpr.NumRows
    for i in 0 .. dim-1 do
        for j in 0 .. dim-1 do
            res <- res + "\t" + DtoStr(bpr.[j,i])
        res <- res + "\n"
    res



(* ------     -----------    --------- --------- --------- --------- *)
(* PR operation definitions                                          *)
(*                                                                   *)
(* We will use the character ^^ to mark PR operators                 *)
(* ------     -----------    --------- --------- --------- --------- *)

let inline (<=^^) (pr1 : PR) (pr2 : PR) : bool = List.forall (fun bpr1 -> List.exists (fun bpr2 -> bpr1 <=^  bpr2) pr2) pr1
let inline (==^^) (pr1 : PR) (pr2 : PR) : bool = (List.length pr1).Equals(List.length pr2) && pr1 <=^^ pr2 && pr2 <=^^ pr1

let rec PRtoSTR(pr : PR) : string = 
    match pr with
        | (bpr::pr_)    -> "BPR:\n" + BPRtoSTR(bpr) + PRtoSTR(pr_)
        | ([]      )    -> "\n"




(* ------     -----------    --------- --------- --------- --------- *)
(* Results pretty-print                                              *)
(* Results pretty-print will extract every variable relation degree  *)
(* from every BPR across all the results to one "Max" degree relation*)
(* Example: 
        PR = 
            BPR1:
	            1	0	0
	            0	1	0
	            0	e	1
            BPR2:
	            1	0	0
	            0	e	0
	            0	1	1
    Variable Xi is represented by the column i of the two BPRs.
    Variable X1 degree relation is the APBs - (1,0,0), (1,0,0) = (1,0,0) = Max(X1) -> X1'~Max(X1).
    Variable X2 degree relation is the APBs - (0,1,e), (0,e,1) = Max(X2, X3) -> X2'~Max(X3, X2).


    Definition:
    Clean AMB - AMB that Epsilons were cleaned, meaning - Epsilon->0
    
    Types:
    CleanAmb   - will represent a clean AMB.
    VarResults - Will contain collection of clean AMBs of single variable Xi from all BPRs.
    Results    - Will contain all Xi VarResults.

    Example - 
            BPR1:
	            1	0	0
	            0	1	0
	            0	e	1
            BPR2:
	            1	0	0
	            0	e	0
	            0	1	1
            -> Results :
                ( ( (1,0,0), (1,0,0) ),     VarResults of X1
                  ( (0,1,e), (0,e,1) ),     VarResults of X2
                  ( (0,0,1), (0,0,1) ) )    VarResults of X3
                    ^CleanAmb1
                             ^CleanAmb2
                                                                     *)
(* ------     -----------    --------- --------- --------- --------- *)
  
type CleanAmb   = int list
type VarResults = CleanAmb list
type Results    = VarResults list

let mutable program : string = ""

let DisclaimerNote = "* Disclaimer:\n    Some of the variables degree relation might be exponential and will only take into consideration at the end of the calculation."

let ResultsPrettyPrint(pr : PR, nonbounded : int list) : string =

    (* Remove Epsilon from D, get a clean D *)
    let cleanD(d : D) : int =
        match d with 
            | (Degree(d)  )   -> int d
            | (Epsilon    )   -> 0

    (* CleanAmb order *)
    let inline GreaterEqualVarResult(vr1 : CleanAmb)(vr2 : CleanAmb) : bool = List.forall2 (>=) vr1 vr2
 
    (* Get rid of unneeded results, for example if x1'=x1 and x1'=x1*x2 then the x1'=x1 is not needed  *)
    let inline MaxResults(results : Results) : Results = List.map (fun vrs -> Max(vrs, GreaterEqualVarResult)) results

    (* Variable and multiplier to string, for example - X1^3 -> X1*X1*X1 *)
    let rec VarMulToStr(varMul : int, varID : int) :string =
        match varMul with
            | 0     -> ""
            | 1     -> "X" + string varID
            | mul   -> "X" + string varID + "*" + VarMulToStr(mul-1,varID)
     
    (* Result to string such as [1,0,3] -> X1*X3*X3*X3 *)   
    let rec VarResultToStr(cleanAmb : CleanAmb, firstVar : int, first_print : bool) : string =
        match cleanAmb with
            | []    -> ""
            | varMul::varResult_   -> (if (first_print || varMul.Equals(0)) then "" else "*") + 
                                        VarMulToStr(varMul, firstVar) + 
                                        VarResultToStr(varResult_, firstVar+1, first_print && varMul.Equals(0))            

    let rec VarResultsToStr(varResults : VarResults) : string =
        match varResults with 
            | ([])                     -> ""
            | (cleanAmb::[])          -> VarResultToStr(cleanAmb, 1, true)
            | (cleanAmb::varResults_) -> VarResultToStr(cleanAmb, 1, true) + ",  " + VarResultsToStr(varResults_)

    (* Check if Variable is nonbounded and should be treated as exponential *)
    let inline IsExp(nonbounded : int list, varId : int) : bool = List.exists varId.Equals nonbounded
    let rec IsExpResult(cleanAmb : CleanAmb, nonbounded : int list, firstVar : int) : bool =
        match cleanAmb with 
            | []    -> false
            | varMul::varResult_ -> ((not (varMul.Equals(0))) && IsExp(nonbounded, firstVar)) || IsExpResult(varResult_, nonbounded, firstVar+1)
    let inline IsExpResults(varResults : VarResults, nonbounded : int list) : bool = List.exists (fun cleanAmb -> IsExpResult(cleanAmb, nonbounded, 1)) varResults
 
    (* Prepare a result, for example :  *)
    (* [[1,0,3],[0,1,0]]                *)
    (*    ->   X1' ~ Max(X1*X3*X3*X3)   *)
    (*         X2' ~ Max(X2)            *)
    let rec ResultsToStr(results : Results) : string =
        let mutable str : string = ""
        for i in 0 .. (results.Length)-1 do
            if (nonbounded.Length > 0 && IsExpResults(results.[i], nonbounded)) then 
                str <- str + "X" + string (i+1) + " is Exponential\n"
            else
                str <- str + "X" + string (i+1) + "' ~ Max(" + VarResultsToStr(results.[i]) + ")\n"
        str
         
    (* Convert RowVector to CleanAmb and by the way cleanD *)
    let inline RowVecToVarResult(rv : RowVector<D>) : CleanAmb =  List.map (fun d -> cleanD(d)) (Array.toList (RowVector.Generic.toArray<D> rv))  


    let numVars = pr.[0].NumCols
    let mutable results = Array.create numVars []
    let mutable res = ""
    for i in 0 .. (List.length pr)-1 do
        for j in 0 .. numVars-1 do
            results.[j] <- RowVecToVarResult(pr.[i].Row j)::results.[j]
    ResultsToStr(MaxResults(Array.toList results))
    
    (* end of ResultsPrettyPrint *)



(* ------     -----------    --------- --------- --------- --------- --------- --------- *)
(* ------     -----------    --------- --------- --------- --------- --------- --------- *)
(* ------     Here is where the Polynomial Complexity algorithm starts         --------- *)
(* ------     -----------    --------- --------- --------- --------- --------- --------- *)
(* ------     -----------    --------- --------- --------- --------- --------- --------- *)

let debugInfoAdded = new Event<DebugInfoAddedEventArgs>()

(* pr_res is actually the algorithm result as PR *)
let mutable pr_res  : PR  = []
let mutable numVars : int = 0
    

type PolyCompAnalysis(analysis_params : Analysis.ANALYSIS_PARMS) =
    let Debug = ref false
    

    (* Event definition for step by step run *)
    [<CLIEvent>]
    member x.CommandGraphUpdated = commandGraphUpdated.Publish

    (* Event definition for debug information *)
    [<CLIEvent>]
    member x.DebugInfoAdded = debugInfoAdded.Publish


    (* Get algorithm results as string using ResultsPrettyPrint *)
    member x.GetResult(nonbounded : int list) : string = 
         "Program:\n\n" + program + "---------------------------\n\nDegree Polynomial Relation =\n\n" + ResultsPrettyPrint(pr_res, nonbounded)

         
    
    (* the entry point for the polynomial complexity execution. Receives a command of the core language as input (in AST form) and 
        creates PR(pr_res) which represents the algorithm analysis results on the input command. *)

    member x.analyse(c : Core.Cmd) : bool =  

        (* debug mode -  when true - each algorithm stage will be printed, along with its output (and sometimes input). *)                
        (* Sends s to the any debug info added event listener. *)
        let debugPrint(s) = if (!Debug) then debugInfoAdded.Trigger(DebugInfoAddedEventArgs(s))
        let debugPrintG(g) = if (!Debug) then debugPrint(CG.prtEdges(g))
        

        let rec RunPolyComp(c : CoreEx.Cmd) = 

                (* Retrieve the number of variables *)
                let rec NumVarsAsgn(e : Core.Expr) : int =   
                    match e with
                        | Core.Var y            -> y
                        | Core.Plus(e1 , e2)    -> max (NumVarsAsgn(e1)) (NumVarsAsgn(e2))
                        | Core.Times(e1 , e2)   -> max (NumVarsAsgn(e1)) (NumVarsAsgn(e2))
                        | _                     -> 0

                let rec NumVars(CoreEx.PolyCompCmd(c,_) as c_: CoreEx.Cmd : CoreEx.Cmd) : int = 
                    match c with
                        | CoreEx.Asgn(id, e)    -> max id (NumVarsAsgn(e))
                        | CoreEx.Seq(c1,c2)     -> max  (NumVars(c1)) (NumVars(c2))
                        | CoreEx.Choice(c1,c2)  -> max  (NumVars(c1)) (NumVars(c2))
                        | CoreEx.Skip           -> 0
                        | CoreEx.Loop(l, c)     -> max l (NumVars(c))
                    

                (* ------     -----------    --------- --------- --------- --------- --------- *)
                (* Expressions to APB conversion, According to inference rules for expressions *)
                (* ------     -----------    --------- --------- --------- --------- --------- *)
                let rec ExpressionsToAPB(e : Core.Expr) : APB =   
                    
                    let zeroList() : AMB = [for i in 1 .. numVars -> Degree(0)]

                    (* Add epsilon to apb1 according to apb2, meaning {amb1 + eps*amb2 | amb1 in apb1 and amb2 in apb2}    *)
                    let rec epsToAPB(apb1 : APB, apb2 : APB) : APB =
                        match apb1 with 
                            | []            -> []
                            | amb1::apb1_   -> List.append (List.map (fun amb2 -> amb1 +& (Epsilon *& amb2)) apb2) (epsToAPB(apb1_, apb2))

                    match e with
                        | Core.Var y                         -> MakeAPBfromAMB(ambSet(zeroList(), y-1, Degree(1)));
                        | Core.Plus(e1 , e2)                 -> 
                            let apb1 = ExpressionsToAPB(e1)
                            let apb2 = ExpressionsToAPB(e2)
                            epsToAPB(apb1, apb2) |=| epsToAPB(apb2, apb1)

                        | Core.Times(e1 , e2)               -> 
                            let apb1 = ExpressionsToAPB(e1)
                            let apb2 = ExpressionsToAPB(e2)
                            apb1 *% apb2                                      
                        | _ -> Util.fail("PolyCompAnalysis.ExpressionsToAPB: Invalid expression.")


                (* --------- ---------    Assignment   --------- --------- --------- --------- --------- *)
                (* Assignment to PR conversion, According to inference rules for commands                *)
                (* Associate PR with every core language command                                         *)
                (*                                                                                       *)
                (* Assignment is evaluating expression into APB and associating the APB with PR          *)
                (* according to assignment variable while using PRset function.                          *)
                (* ------     -----------    --------- --------- --------- --------- --------- --------- *)
                let Asgn(e : Core.Expr, id : int) : PR = 
                    (* PRset is APB to PR conversion that associates every AMB with BPR.           *)
                    (* For every AMB of APB - Set the AMB at the #id column of the new BPR, Every  *)
                    (* other column are identity in BPR matrix.                                    *)
                    let rec PRset(apb : APB, id : int) : PR =
                        let inline BPRset(amb : AMB, id : int) : BPR = 
                            Matrix.Generic.mapi<D> (fun i j d -> if i.Equals(id-1) then amb.[j] else d) (IdentityMatrix(numVars))
                        List.map (fun amb -> BPRset(amb, id)) apb    
                    
                    PRset(ExpressionsToAPB(e), id)


                (* --------- --------- Choice --------- --------- --------- *)
                (* PRs Union, According to inference rules of commands      *)
                let Choice(pr1 : PR, pr2 :PR) : PR =
                    Max(List.append pr1 pr2, (>=^))


                (* --------- --------- Sequence --------- --------- --------- *)
                (* Multiply PRs, According to inference rules of commands     *)
                (* Multiply every bpr1 in pr1 with every bpr from pr2.        *)
                let rec Seq(pr1 : PR, pr2 :PR) : PR =
                    let inline SeqBPR(bpr : BPR, pr :PR) : PR = List.map (fun bpr_ -> bpr *^ bpr_) pr
                    match pr1 with 
                        | []           ->  []
                        | bpr1::pr1_   ->  List.append (SeqBPR(bpr1, pr2)) (Seq(pr1_, pr2))  


                

                (* --------- --------- Loop --------- --------- *)
                (* The process of the loop consists of the closure, loop correction and again closure. *)
                (* The closure should consider exponential variables that might be exist in loop command scope in order to avoid infinite loop of closure *)

                (* GetExpFlags returns a Boolean list of variable flags while flags[i]=true if and only if Xi is exponential variable if PR. *)
                (* The columns of exponential will be ignored in Closure *)
                (* Xi is direct exponential if 
                       BPR[i,i]>1 such as Xi'=Xi*Xi
                    or BPR[i,i]==1 and exist j!=i such that BPR[i,j]>=1 such as Xi'=Xi*Xj 
                   Xi is indirect exponential if 
                    Xi depend on Xj that is exponential such as Xj'=Xj*Xj and Xi'=Xj *)
                let GetExpFlags(pr : PR) : bool list = 
                    (* Find direct exponential, in every function call we get one flag of current column(by the argument 'var') and recursively get the rest of the list *)
                    let rec GetDirectExp(bpr : BPR, var : int) : bool list = 
                        let varIsExp = 
                            ( bpr.[var,var] >?  Degree(1) ||       (* BPR[i,i]>1 ? *)
                              (bpr.[var,var] ==? Degree(1) && Vector.Generic.existsi<D> (fun j d -> j<>var && d >=? Degree(1)) (bpr.Row var).Transpose )  )  (* exist j such that BPR[i,j]>=1 ? *)
                        let restOfTheList = (if (var < numVars-1) then GetDirectExp(bpr, var+1) else []) 
                        (varIsExp :: restOfTheList)
                
                    (* Find indirect exponential *)
                    let rec GetIndirectExp(bpr : BPR, expFlags : bool list) : bool list =
                        (* In every function call we get recheck one flag by current column(be the argument 'var') and recursively check the rest of the list *)
                        let rec CalcIndirectExp(bpr : BPR, directExp : bool list, var : int) : bool list = 
                            let varIsExp = 
                                ( directExp.[var] ||         (* Direct exponential *)
                                  Vector.Generic.existsi<D> (fun j d -> d >=? Degree(1) && directExp.[j]) (bpr.Row var).Transpose) (* Indirect exponential *)
                            let restOfTheList = (if (var < numVars-1) then CalcIndirectExp(bpr, directExp, var+1) else []) 
                            List.append [varIsExp] restOfTheList

                        (* Continue searching while finding new indirect exponential variables *)
                        let next_expFlags = CalcIndirectExp(bpr, expFlags, 0)
                        if (expFlags <> next_expFlags) then
                            GetIndirectExp(bpr, next_expFlags)
                        else
                            expFlags

                    (* Get exponential variables list *)
                    let GetExp_BPR(bpr : BPR) : bool list = 
                        GetIndirectExp(bpr, GetDirectExp(bpr, 0));

                    (* Perform or operation between all exponential variables lists according to list index, example: [(f,t,f),(f,f,f),(f,f,t)]->(f,t,t) *)
                    let rec ListOR(ll : bool list list, var : int) : bool list = 
                        ( List.exists (fun l -> (List.nth l var)) ll ) :: (if (var < numVars-1) then ListOR(ll, var+1) else []) 

                    ListOR( (List.map GetExp_BPR pr), 0)

                    (* end of GetExpFlags *)

                (* Check if pr1<=pr2 while ignoring rows of exponential variables *)
                let SmallerOrEqual_PR_IgnorExp(pr1 : PR, pr2 : PR, expFlags : bool list) : bool =
                    (* bpr1<=bpr2 if all D in bpr1 are smaller/equal matching D in bpr2 or D is in exponential variable row *)
                    let inline SmallerEqual_BPR_IgnorExp      (bpr1 : BPR, bpr2 : BPR, expFlags : bool list) : bool = Matrix.Generic.foralli (fun i j d1 -> d1 <=? bpr2.[i,j] || expFlags.[i] ) bpr1
                    let inline SmallerEqual_BPRfromPR_IgnorExp(bpr1 : BPR, pr2  : PR , expFlags : bool list) : bool = List.exists (fun bpr2 -> SmallerEqual_BPR_IgnorExp(bpr1, bpr2, expFlags)) pr2
                    List.forall (fun bpr1 -> SmallerEqual_BPRfromPR_IgnorExp(bpr1, pr2, expFlags)) pr1
              
                (* Closure is the smallest set that contains pr and closed under multiplication *)
                (* Closure will take into consideration exponential variables while ignore their columns *)
                let rec Closure(pr : PR) : PR =

                    let () = debugPrint("--------------------------------\n\nPR for Closure = \n");
                    let () = debugPrint(PRtoSTR(pr));
                
                    let expFlags = GetExpFlags(pr);
                    let pr_next = List.append pr (Seq(pr, pr));
                    (* Check only if pr_next<=pr because pr is part of pr_next[pr_next>=pr always true] *)
                    if (SmallerOrEqual_PR_IgnorExp(pr_next, pr, expFlags)) then
                        pr 
                    else 
                        Closure(pr_next)

                (* Correct a PR at variable l according to idempotent and epsilon on the main diagonal *)
                let LoopCorrection(l : int, pr : PR) : PR =

                    let () = debugPrint("--------------------------------\n\nPR for LoopCorrection = \n");
                    let () = debugPrint(PRtoSTR(pr));
                    let () = debugPrint("--------------------------------\n\n");

                    (* if bpr[i,i]=epsilon then we need to correct Xi by multiply Xi by Xl *)
                    let CorrectD_in_BPR = fun (bpr : BPR) (i : int) (j : int) (d : D) (l : int) ->  
                        if (bpr.[i, i].Equals(Epsilon)) then    (* Correct the column only the variable at the main diagonal is epsilon *)
                            if (j.Equals(l-1)) then
                                d +? Degree(1)      (* Correction! by multiplying by Xl *)
                            else 
                                d   (* no correction is needed *)
                        else 
                            d   (* no correction is needed *)


                    (* Go over all BPRs of PR and correct the idempotent BPRs *)
                    let CorrectBPRifIidempotent = fun (bpr : BPR) (l : int) (pr_num : int) ->      
                        (* Go over all entries in BPR and check for correction *)
                        let CorrectBPR(l : int, bpr : BPR) : BPR = Matrix.Generic.mapi<D> (fun i j d -> (CorrectD_in_BPR bpr i j d l) ) bpr
                    
                        (* Is idempotent? (BPR*BPR==BPR)  *)
                        if (Seq([bpr], [bpr]) ==^^ [bpr]) then (* Correct the BPR only if idempotent matrix *)
                            let pr_corrected = CorrectBPR(l, bpr)      
                            let () = debugPrint("Corrected PR[" + string (pr_num+1) + "] = \n" + BPRtoSTR(pr_corrected));
                            pr_corrected
                        else
                            let () = debugPrint("PR[" + string (pr_num+1) + "] = is not idempotent, Correction is not needed\n");
                            bpr (* no correction is needed *)
                    
                    (* Correct all BPRs in pr *)
                    List.mapi (fun i bpr -> (CorrectBPRifIidempotent bpr l i) ) pr
                    
                    (* end of LoopCorrection *)

                let Loop(l : int, pr : PR) : PR =
                    Max(Closure(LoopCorrection(l, Closure(IdentityMatrix(numVars)::pr))), (>=^))

                    




                (* --------- ---------  --------- --------- --------- --------- --------- --------- *)    
                (* --------- --------- Recursive PR evaluate from core language --------- --------- *)    
                (* --------- ---------  --------- --------- --------- --------- --------- --------- *)
                (* This is the algorithm main function that recursively analyze the program and     *)
                (* returning the pr that represent the program                                      *)
                (* --------- ---------  --------- --------- --------- --------- --------- --------- *)
                let rec PREvaluate(CoreEx.PolyCompCmd(c,_) as c_: CoreEx.Cmd : CoreEx.Cmd) : PR = 
                    let pr = match c with
                            | CoreEx.Asgn(id, e)    -> Asgn(e, id)
                            | CoreEx.Seq(c1,c2)     -> Seq(PREvaluate(c1),PREvaluate(c2))
                            | CoreEx.Choice(c1,c2)  -> Choice(PREvaluate(c1),PREvaluate(c2))
                            | CoreEx.Skip           -> [IdentityMatrix(numVars)]
                            | CoreEx.Loop(l, c)     -> Loop(l, PREvaluate(c));
                    (* send event that command graph was updated *)
                    commandGraphUpdated.Trigger(CommandGraphUpdatedEventArgs(c_, null, "Degree Polynomial Relation =\n\n" + ResultsPrettyPrint(pr, []) + "\n" + DisclaimerNote + "\n\n---------------------------\n\nPR =\n\n" + PRtoSTR(pr)))
                    pr
                    
                    
                (* --------- ---------  --------- --------- --------- --------- *)    
                (* --------- --------- Run the algorithm    --------- --------- *)    
                (* --------- ---------  --------- --------- --------- --------- *)    
                let () = debugPrint("starting to analyze command:\n");
                program <- CoreEx.prt(c)
                let () = debugPrint(program);
                (* Compute the number of variables *)
                numVars <- NumVars(c);
                
                (* Run algorithm and save to pr_res *)
                pr_res <- PREvaluate(c);
                
                let () = debugPrint("Finished analyzing program, result PR = \n");
                let () = debugPrint(PRtoSTR(pr_res));
                true


        (* --------- ---------  --------- --------- --------- --------- --------- *)    
        (* --------- --------- Entry point of algorithm execution       --------- *)    
        (* --------- ---------  --------- --------- --------- --------- --------- *)    
        RunPolyComp(CoreEx.extendedPolyCompAST(c))
    
        (* end of x.analyse *)
               

    (* set the debug output mode *)
    member x.setDebugMode(mode : bool) = Debug := mode
    
    (* end of PolyCompAnalysis *)


let Poly = PolyCompAnalysis(Analysis.PolyParams())      